Context

One of the biggest causes of mortality among women is breast cancer [1]. The more the detection is delayed, the lower the chance of survival is. That is why efficient diagnosis is critical to provide personalized care and treatment, and to improve cancer outcomes. However, it is quite difficult due to the complexity and diversity of this disease.

A major concern of the field is to be able to classify the type of tumor with very good accuracy and determine critical features that allow this classification. Particular attention is paid to the false positive rate of this classification since it is a vital issue.

1. Information about the dataset

We propose to study the Wisconsin Breast Cancer Dataset.

To understand the feature of the selected dataset, we will need first some definitions:

• According to [2], a tumor is “An abnormal mass of tissue that forms when cells grow and divide more than they should or do not die when they should. Tumors may be benign (not cancer) or malignant (cancer). Benign tumors may grow large but do not spread into, or invade, nearby tissues or other parts of the body. Malignant tumors can spread into, or invade, nearby tissues”

• Specifically, breast cancer is a malignant tumor located in breast cells of the fatty tissues or the fibrous connective tissues within the breast. [3]

• “Fine needle aspiration is a type of biopsy procedure. In fine needle aspiration, a thin needle is inserted into an area of abnormal-appearing tissue or body fluid.” [4]

Benign and malignant tumor of a breast mass, fine needle aspiration [5]

According to data description, “the features are computed from a digitized image of a fine needle aspirate (FNA) of a breast mass. They describe characteristics of the cell nuclei present in the image.”

There are ten real-valued features computed for each cell nucleus:

The mean, standard error, and “worst” or largest (mean of the three largest values) of these features were computed for each image, resulting in 30 features. All feature values are recorded with four significant digits.

Note: The following analysis has been done with the full dataset. However, the mean of each feature has shown to be more relevant, so we decide to present the analysis only with the mean variables.

Link to dataset

2. Question of interest

Find a model (classifier) able to predict whether the tumor is benign or malignant with only a few features.

3. Plan of the analysis

• Inspect the data to understand it better
• Check if the sample proportion of benign and malignant cancer is representative of the whole population
• Select variables of interest by discarding according to the correlation level
• Build GLM models and select the best model with anova, question the validity of this model
• Address the question of accuracy of the prediction and false positive rate of the selected model
• Apply PCA to see to what extent dimension can be reduced and perform a GLM model with relevant dimensions.

In the annex, we build a GLM model with the library “glmnet” and we find a similar selection of variables.

1. Required libraries

# Check if packages are installed and if not install them
if(!require(pacman)) {
  install.packages(c("pacman", "remotes", "equatiomatic"))
}
if (!require(papaja)) {
  remotes::install_github("crsh/papaja")
  remotes::install_github("kassambara/factoextra")
}
# if(!require(devtools)) {
#   install.packages("devtools")
#   
# }
# if(!require(caret)) {
#   install.packages("caret")
# }
suppressPackageStartupMessages(pacman::p_load(pander,       # Rmarkdown visualization
               GGally, 
               ggfortify,
               ggplot2,
               here,         # to load the path
               kableExtra,   # Rmarkdown table visualization
               papaja,
               glmnet,       # GLM implementation and analysis library
               equatiomatic, # Rmarkdown model equation visualization
               patchwork,    # arrange subplot
               devtools,     # tool for PCA
               factoextra,   # tool for PCA
               caret         # for the confusion matrix
               
))
knitr::opts_chunk$set( comment = NA)

2. Import data

path = here("LUCILE") # get relative path
setwd(path) # set working directory
df <- read.csv('data.csv', stringsAsFactors = 1)# load data 

Let’s delete ID number and the last variable - which is full of NA- because there are not relevant.

df<-df[,-33]; df<-df[,-1]

We will work only with the mean-type data, so let’s create a new frame for the variables of this type.

df_mean <- data.frame(
  "diagnosis"          = df$diagnosis,
  "radius"             = df$radius_mean,
  "texture"            = df$texture_mean,
  "perimeter"          = df$perimeter_mean,
  "area"               = df$area_mean,
  "smoothness"         = df$smoothness_mean,
  "compactness"        = df$compactness_mean,
  "concavity"          = df$concavity_mean,
  "concave.points"     = df$concave.points_mean,
  "symmetry"           = df$symmetry_mean,
  "fractal_dimension"  = df$fractal_dimension_mean
)

4. Proportion of benign vs malignant cancer

kable(prop.table(table(df_mean$diagnosis)),col.names = 
        c("Type","Frequency"))   

The two types of tumors are not represented in the same proportion, this can lead to a bias.

However, according to [6] , this proportion seems to be more or less representative of the reality:

“The benign to malignant ratio (B:M ratio) among breast biopsies (number of benign breast lesions divided by number of breast cancers) is widely believed to be around 4:1 or 5:1”

2. Variables selection

• As expected, radius, perimeter and area are highly correlated (\(r \approx 1\))

• Surprisingly, concavity, compactness and concave.points have a strong correlation.( \(r \approx 0.8\) or \(r \approx 0.9\))

Note:

Even if compactness is defined as \(perimeter^2 / area - 1.0\) the r between this variable and area or perimeter is not 1 because the correlation shows only the linear dependency and their relationship is not linear.

m_perimeter <-
  lm(
    data = df_mean,
    perimeter ~ radius + texture + area + smoothness + compactness + concavity +
      concave.points + symmetry + fractal_dimension
  )
pander(summary(m_perimeter))

m_area <-
  lm(
    data = df_mean,
    area ~ radius + texture + perimeter + smoothness + compactness + concavity +
      concave.points + symmetry + fractal_dimension
  )
pander(summary(m_area))

m_compactness <-
  lm(
    data = df_mean,
    compactness ~ radius + texture + perimeter + smoothness + area + concavity +
      concave.points + symmetry + fractal_dimension
  )
pander(summary(m_compactness))

df_mean_reduc <-df_mean[-c(4,5,7)]

1. Set a first GLM model

We set a GLM model with the remaining features. Since we want to classify two types of tumors, we use the binomial family.

m <-
  glm(
    data = df_mean_reduc,
    diagnosis ~ radius + texture + smoothness + concavity + concave.points +
      symmetry + fractal_dimension,
    family = binomial
  )
summary(m)

Call:
glm(formula = diagnosis ~ radius + texture + smoothness + concavity + 
    concave.points + symmetry + fractal_dimension, family = binomial, 
    data = df_mean_reduc)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.35180  -0.13938  -0.03229   0.02046   3.15368  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -28.38387    6.66946  -4.256 2.08e-05 ***
radius               0.88701    0.21852   4.059 4.93e-05 ***
texture              0.37262    0.06212   5.998 2.00e-09 ***
smoothness          78.50170   32.64920   2.404   0.0162 *  
concavity           15.52082    8.35462   1.858   0.0632 .  
concave.points      46.67203   26.16265   1.784   0.0744 .  
symmetry            16.85783   10.75613   1.567   0.1170    
fractal_dimension -101.54448   61.26233  -1.658   0.0974 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 751.44  on 568  degrees of freedom
Residual deviance: 153.35  on 561  degrees of freedom
AIC: 169.35

Number of Fisher Scoring iterations: 8

• The algorithm converges: the number of fisher scoring iterations is reasonable.
• There is a lot of variables not significant. To solve this we can try to remove them.
• The ratio of the residual deviance by its degrees of freedom is \(153/561 = 0.272\) where the dispersion parameter is 1. There is underdispersion. To solve this we can use the quasibinomial family.

3. Deal with underdispersion

Let’s set a GLM model with a quasibinomial family to solve the issue of underdispersion.

m_selected_quas <-
  glm(
    data = df_mean_reduc,
    diagnosis ~ radius + texture + smoothness + concave.points,
    family = quasibinomial
  )
(summary(m_selected_quas))

Call:
glm(formula = diagnosis ~ radius + texture + smoothness + concave.points, 
    family = quasibinomial, data = df_mean_reduc)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.42132  -0.15010  -0.04247   0.02603   2.86598  

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -28.57552    3.25917  -8.768  < 2e-16 ***
radius           0.85081    0.11585   7.344 7.30e-13 ***
texture          0.35845    0.04052   8.847  < 2e-16 ***
smoothness      52.26403   17.65979   2.959  0.00321 ** 
concave.points  78.73692   11.23385   7.009 6.88e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasibinomial family taken to be 0.458343)

    Null deviance: 751.44  on 568  degrees of freedom
Residual deviance: 160.32  on 564  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 8

The formula of the retained model is:

extract_eq(m_selected_quas, use_coefs = TRUE, wrap = TRUE, terms_per_line = 2)

\[ \begin{aligned} \log\left[ \frac { \widehat{E( \operatorname{diagnosis} )} }{ 1 - \widehat{E( \operatorname{diagnosis} )} } \right] &= -28.58 + 0.85(\operatorname{radius})\ + \\ &\quad 0.36(\operatorname{texture}) + 52.26(\operatorname{smoothness})\ + \\ &\quad 78.74(\operatorname{concave.points}) \end{aligned} \]

The ratio of the residual deviance by its degrees of freedom is \(160.32/564 = 0.284\) where the dispersion parameter is 0.458. These two values are close, we can validate this model.

See in annex another GLM modelization with glmnet.

4. Confusion matrix

First, let’s separate the data into train and test sets.

# set the seed to make partition reproducible
set.seed(6)
prop_train_test <- floor(0.75 * nrow(df_mean_reduc))

# 0.75 of the data fall randomly into the training set and the remaining is for the test set
train_ind <- sample(seq_len(nrow(df_mean_reduc)), size = prop_train_test) 
train <- df_mean_reduc[train_ind, ]
test <- df_mean_reduc[-train_ind, ]
y_test<- as.numeric(test$diagnosis)-1

Built the model with the data of the train set

m_selected_quas_test <-glm(data = train,
                          diagnosis ~ 
                          radius + texture + smoothness + concave.points,
                          family = quasibinomial)

Test the accuracy of the model with the test set by computing the confusion matrix and other statistics.

pdata <- predict(m_selected_quas_test, newdata = test, type = "response")
pred <- as.numeric(pdata>0.5)
conf_mat<-confusionMatrix(reference = as.factor(y_test), data = as.factor(pred))
rownames(conf_mat$table) =c("B","M")
colnames(conf_mat$table) =c("B","M")
kable(prop.table(conf_mat$table))

• True positive and false negative are high (\(TP \approx 0.59\),\(FN \approx 0.36\))
• True negative and false positive are low ($ TN $, $ FP $ )

conf_mat

Confusion Matrix and Statistics

          Reference
Prediction  B  M
         B 84  4
         M  4 51
                                          
               Accuracy : 0.9441          
                 95% CI : (0.8927, 0.9755)
    No Information Rate : 0.6154          
    P-Value [Acc > NIR] : <2e-16          
                                          
                  Kappa : 0.8818          
                                          
 Mcnemar's Test P-Value : 1               
                                          
            Sensitivity : 0.9545          
            Specificity : 0.9273          
         Pos Pred Value : 0.9545          
         Neg Pred Value : 0.9273          
             Prevalence : 0.6154          
         Detection Rate : 0.5874          
   Detection Prevalence : 0.6154          
      Balanced Accuracy : 0.9409          
                                          
       'Positive' Class : 0               
                                          

The accuracy of the model is 0.9441, with a 95% confidence interval:(0.8927, 0.9755).

Further direction: Here we compute the confusion matrix and the other statistics with only one test set. It could be accurate to do this computation for several test set.

2. Observations in PC plans

autoplot(p,x = 1,y = 2, data = df_mean_reduc,colour = "diagnosis",
  loadings = TRUE, loadings.colour = "blue", loadings.label = TRUE)+
  theme_minimal()

By plotting the data in the plan of the first 2 principal components we see a clear separation between benign and malignant type of cancer. That mean that knowing the location in this plan of a new observation should allow to predict cancer. Let’s try this approach by applying a GLM model to the first two PC.

Note:

We tried to plot the data in all the possible pairs formed by the first four components. The clearest separation between the type of cancer occurs in the plan PC1-PC2.

3. GLM model with PCA

df_mean_pca <- cbind(df_mean_reduc, p$x)
glm_pca <-glm(data= df_mean_pca, df_mean_pca$diagnosis ~PC1+PC2, 
              family = quasibinomial)

summary(glm_pca)

Call:
glm(formula = df_mean_pca$diagnosis ~ PC1 + PC2, family = quasibinomial, 
    data = df_mean_pca)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.74789  -0.16378  -0.04399   0.02463   3.13937  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.6441     0.1600  -4.026 6.44e-05 ***
PC1          -3.1546     0.2727 -11.569  < 2e-16 ***
PC2           2.4026     0.2384  10.077  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasibinomial family taken to be 0.6184117)

    Null deviance: 751.44  on 568  degrees of freedom
Residual deviance: 162.50  on 566  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 8

1.diagnostic plot of retained glm model

 autoplot(m_selected_quas,1:4)+ theme_minimal()

* The “residual vs fitted” and “scale-location” subplots show a local regression curve close to 0.

• In the “Q-Q plot,” the central points seem to follow well the normality but it is not the case at the extremes of the distribution.

• In the “Cook’s distance” plot there is no observation with a high influence on the model: all distances are smaller than 0.5.

These observations allow us to conclude that this model can be validated.